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Abstract 

We study a two-type branching process which provides excellent description of experimental data on cell dynamics 
in skin tissue (Clayton et al. |2007[ l. The model involves only a single type of progenitor cell, and does not require 
support from a self-renewed population of stem cells. The progenitor cells divide and may differentiate into post- 
mitotic cells. We derive an exact solution of this model in terms of generating functions for the total number of cells, 
and for the number of cells of different types. We also deduce large time asymptotic behaviors drawing on our exact 
results, and on an independent diffusion approximation. 



1. Introduction 



Understanding the kinetics (homeostasis) of cells in adult mammalian tissues has long been a major challenge in 
biology. Recent progress in experimental methods made it feasible to label individual cells in vivo, and follow their 



fate and that of their progeny (Clarke and Tickle 1999 Jonkers and Berns 2002 1. This powerful genetic labeling 
technique has enabled in vivo experiments in the outmost layer of skin (epidermis) of the tail in adult mice ( !Clayton| 
|eral.,,2007^) . Individual cells in the basal layer of the epidermis have been marked by a fluorescent genetic label and 
the size of the clone (all living progenies of a cell) of each single marked cell has been measured at different times. 
This has provided the data on the evolution of the clone size distribution in the basal layer of the epidermis. 

The prevailing model of epidermal homeostasis has involved long-lived stem cells generating short-lived popu- 
lations of transit-amplifying (TA) cells that differentiate into post-mitotic cells ( |Potten[ |1974[ [Blanpain and Fuchs} 
2009| l. The stem-TA hypothesis predicts that the clones of TA cells should disappear (after sufficiently long time). 



while the existing clones should be small and associated with stem cells. Strikingly, the fraction of remaining clones 
was found to decrease as (time) accordingly the average size of existing clones scales linearly with time. This re- 



markable scaling behavior calls for a totally different model of epidermal homeostasis. Clayton et al. (2007 1 proposed 
a model of cell division and differentiation which manifestly obeys the observed scaling behavior and provides an 
excellent fit to more subtle characteristics. A gratifying property of the model suggested by [Clayton et al.| ( |2007| l is 
that it is simpler than the stem-TA model: The new model involves only a single type of committed progenitor cell, 
and in particular, stem-cell proliferation is not required for epidermal homeostasis. 

Thus the model describes the population of cells of two types: Proliferating cells (type A) divide and eventually 
differentiate into non-proliferiting cells (type B), which leave the basal layer and migrate to the epidermal surface 
where they are shed. More precisely, the cell population evolves according to the continuous time, constant rate, 
two-type branching process 



A 
A 
A 
B 



AA 
AB 
BB 

■0 



at rate r 
at rate \ -2r 
at rate r 
at rate y 



(1) 
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Here we set the overall cell division rate to unity. In the experiments of Clayton et al. ( 2007| l, the division rate was 
equal to A - 1.1 /week; the values of the parameters were found to be r = 0.08 and y = F/A = 0.28. Note that the 
model is assumed to be critical, that is the division rates corresponding to the channels A — > AA and A — > BB are the 
same. Due to this symmetry, the average population size of progenitor cells remains constant as it is required by the 
steady-state assumption. The average population size of the post-mitotic cells is also constant. 

So far the model has been experimentally tested only in mice tail skin. There are still technical constraints pre- 
venting the quantitative tests of the model in other tissues, but those problems are temporary. The model challenges 
the necessity of stem-cell proUferation for the homeostasis of epidermis ( [Jones and Simons] |200 8 ). There is also 
growing evidence ( Dor et al. 2004 Giangreco et al. 2009^ that stem cells do not contribute to the maintenance of 
various other adult tissues. Hence the two-type branching process ([T]) may find a broad range of applications and 
therefore it is highly desirable to possess an exact solution. Despite its apparent simplicity, the branching process 
([T]i has not been solved, although some exact and asymptotic behaviors have been found ( Clayton et al. 2007 Klein 
|et al.|[2007l|2008| l. In this paper we apply generating function techniques to obtain an exact analytic solution, as well 
as approximate methods to derive asymptotic limits. 

Branching processes have been extensively used to model proliferation of differentiating cells, especially in the 
hemopoietic (blood production) system ( Vogel et al. 1969 Pharr et al. 1985[ l; see also other references in Section 

12002). An interesting multi-type model has also been proposed recently in (Dingli 



6.9.1 in (Kimmel and Axcelrod 



et al. 2007b,, 2009] l. These studies, however, mainly rely on numerical solutions, while analytic treatment is restricted 



to obtaining average quantities (or second moments). 

The rest of the paper is organized as follows. We introduce the model in Section[2] and discuss its basic behavior. 
After presenting the generating function methods in Section[3] we provide an elementary solution on a special line in 
the parameter space in Section |4] The model admits a neat explicit solution at the special point y - l,r - 1/4, which 
is discussed in Section |5] As our main result, we derive the generating function of the model for general parameter 
values in Section |6] where we also present an efficient numerical method to obtain the probabilities of having certain 
number of cells at a given time. We discuss the large time asymptotic behavior in Section [7] and derive additional 
scaling properties by means of the Fokker-Plank method in SectionjH] Final remarks are presented in Section|9] 



2. The Model 



The model involves two types of cells, A and B. Type A cells (progenitor cells) are able to divide (proliferate) and 
diffirentiate, B cells (post-mitotic cells) do not divide, they just die (leave the basal layer). More precisely, the two 
cell populations evolve according to the two-type branching process ([T]). The probabiUty Pm,n{f) of having m copies 
of A, and n copies of B at time t satisfies 



dt 



= r(m - l)P,„-iM + (1 - 2r)mP,„,„_i + r(m + l)Pm+\,n-2 + y{n + l)^'m,„+i -im + yn)P„,,„ 



(2) 



The consecutive gain terms on the right-hand side of Eq. (|2]i merely describe the contributions of the channels (from 
top to bottom) of the two-type branching process ([T}. To determine the clone size distribution we start with a single A 
cell, that is 

Pm,n(t = 0) = 5,n,ldn,0 (3) 

We are interested in the full distribution Pm,n(t) and also in the reduced probability distribution Tls(t) of having s - m+n 
total cells at time f; the latter distribution is directly probed in experiments. Needless to say. 



(4) 



Let us first determine the probability distribution P,„(0 of having m cells of type A. This probability distribution 
is readily found since B cells do not affect A cells, and A cells alone evolve according to the critical branching process 



A^AA 
A ^ 



at rate r 
at rate r 



(5) 



2 



The solution, for the initial condition P,„{t = 0) = is ( Athreya and Ney 
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Pm(t) = 

1 + rt 

Notice that the average number of A cells remains constant 



1 / rt V"-^ ^ 

^ for m > 1 

+ rtV \ l + rt/ 



(l+rty\l+rtJ (6) 
rt ^ 
for m - 



<m)= ymP,„(f) = 1, (7) 



m>0 

throughout the evolution. This is of course a general property of the critical branching process. 

We can also compute the average number of post-mitotic cells (n) = Yjm,n>o nPm,n{t)- Indeed, this quantity satisfies 
a simple rate equation 

din) 

= <m> - y{n) (8) 

at 

The gain term on the right-hand side of ([8| follows from the second and third channels (from top to bottom) of the 
two-type branching process ([T]l; the loss term corresponds to the last channel. Using (m) = 1 and («)|r=o = we solve 
^ to yield 

1 - e-^' 

in) = (9) 

r 

Therefore the total average number of cells is given by 

1 - e-^' 

{s)^\ + (10) 

r 

Note that the fraction of type A cells is asymptotically 

^^7T = T^ (11) 

{s) 1 + y 

These exact expressions for the average population sizes are useful and e.g. the fraction of type A cells ( [TT| i will 
appear in numerous latter formulae. The full description of the clone size requires analyzing an infinite set of master 
equations (|2|. We shall perform such analysis using generating function techniques. 



3. Generating function 

We define the generating function of Pm,n{f) as 

oo 

F{x,y,t)^ ^ x'«/P„,„(f) (12) 

Note that (|2| is valid for all m,n > 0, if we define Pm^„ = for all m < or n < 0. [Such systems are said to have 
natural boundary conditions (van Kampen 1997| l.] We multiply both sides of (|2| by x^"y" and sum over all values of 
m,n > 0. Using identities mx'" - xdxxf",ny" - ydyy", where = d/dx,dy = d/dy, we arrive at a partial differential 
equation 

d,F +[x(l-y)- r(x - yf] d,F + y{y-l)dyF ^0 (13) 
The initial condition ([3]l corresponding to a single initial A cell becomes 

F(x,y,t^O)^ X (14) 

Thus we need to solve the partial differential equation ( [T3| ) subject to ( [l4j i. Mathematically, equation ( [T3| ) is a 
hyperbolic partial differential equation and it can be analyzed using the method of characteristics (Logan] 2008| l. 
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Instead, we employ backward Kolmogorov equations; this approach is technically somewhat easier in the present 
case. Here we need two generating functions Fa and Fb, where the subscripts refer to the type of the single initial 
cell. For the forward case we only needed the interesting F = Fa- The initial conditions are 



FA(,x,y, t - 0) - X 

(15) 

FB(x,y,t = 0)=y 
The coupled backward Kolmogorov equations read 

d,FA = rFl + (1 - 2r)FAFB + rFl - Fa (16a) 
d,FB^y{\-FB) (16b) 

These equations can be derived from the corresponding backward Kolmogorov equations for the probabilities Pm^„{t), 



or they can be written down directly ( Athreya and Ney 2004j). The negative terms (-Fa, -jFb) describe the disap- 



pearance of a cell, and the positive terms stand for the created new cells, with the corresponding rates. The term "1' 
in \ \6h\ is just the generating function of no created particle, that is 1 = x°y°. 
Equation ( |16b| i is immediately solved to give 

Fb = ye-^' + 1 - e-^' = f (17) 

This is not surprising, of course: Starting with a single B cell, the system will either contain the initial B cell (this 
occurs with probability e ''') or no cells at all. Substituting ( [17] ) into ( |16a| i and changing the variable from f to / we 
obtain 

y(l-f)dfF^r(F-ff+F(f-l) (18) 

where we dropped the subscript A so that F = Fa- We further simplify the above equation by changing variable / to 
u - 1 - f - (1 - y)e^'^' - The function F{u) then satisfies 

yuF' ^ Fu - r(F + u - if (19) 

with initial condition F{u = I - y) - x- In equation ([T9| and later the prime denotes the derivative with respect to m. 
Note that the forward equation ( [T3] l leads to the same equation ( [T9] l via the method of characteristics. 

Equation ([T9| is an ordinary differential equation of the first order. Yet it is non-linear and could be unsolvable 
as it belongs to the family of Riccati equations. Riccati equations are in principle intractable, yet there are two tricks 
which sometimes allow to solve certain Riccati equations (Bender and Orszag 1978|. One is based on the reduction 



to the linear ordinary differential equation of the second order, the Sturm-Li ouville equation. Another trick applies 
if we manage to find a special solution. We shall see that both tricks lead to success. Let us begin with the more 
elementary second approach. 



4. Elementary Solutions 

The idea is to guess one solution Ftiu) irrespective whether it satisfies the initial condition or not. Having found 
such a special solution, one then seeks a general solution in the form 

F{u) = FAu) + (20) 
V(u) 

The function V(u) satisfies a linear differential equation which is readily solvable. 
The form of ^T9[ suggests to seek a special solution as a polynomial: 

Ft(u) = Aq+Aiu + . - - +Apu'' (21) 

Here Aq, Ai, . - - ,Ap are constants and Ap + 0, so that the polynomial ( |2T| ) has degree p- Noting that uF' is the 
polynomial of degree p, uF is the polynomial of degree p + I, and {F - \ + u)^ is the polynomial of degree 2p, 
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equating the highest degree in u would be possible only if /? + 1 - 2p, i.e. p - I - Thus the polynomial solution should 
be a linear function of u, 

F^{u}^Aq+Axu (22) 

Plugging ( |22| ) into ([T9| we find that the matching is achieved [that is, the ansatz ( |22) i works] if Aq = l,Ai = l/y and 
the parameters r, y are related via 

r 



7 



(i+r)2 

The prescription ( |20l ) tells us to seek the general solution in the form 

u 1 

F = 1 + 



7 V(u) 

By inserting ( |24] i into ([T9]l we arrive at a linear ordinary differential equation 



V +yV- 



(1 M 



- = 0, 



. _ 1 r- 1 

^ y y + 1 



The homogeneous part has solution e ^" and therefore the general solution to 
auxiliary function W obeys 

1 e^" 



W 



(1 + 7)2 M 



which is solved to yield 



W = -r + const 

(1 +r)2 



(23) 
(24) 

(25) 

is sought as V = e"^" W. The 

(26) 

(27) 



/-OO J- 

Here EiXjc) = - J d^e^^ is the exponential integral. The constant (and the choice of the appropriate low limit in 
the integral) in ( |27| l are fixed by the initial condition. Recall that initially we have F{u - 1 - y) = x. Hence ( |24] l gives 



1 + 



y-l 



1 

Vo 



and therefore 



Using (|27]l and (|29]l we obtain 



f(i-y) 



where we used the shorthand notation 

&(y,t)-- 

Combining (|24]i and (|30| we arrive at 



e 



y(i-y) 



.X - 1 +(y- l)/7 



.X - 1 + Cy- i)/r 

£'[y(i-y)]-£/[y(i-y)g'n 
(i+r)2 



F(jc,y,f) = 1 +r"'e"'"(l -y) + 



pf(l-y) 



-£(y,f) 



(28) 



(29) 



(30) 



(31) 



(32) 



x-l+iy- l)/y 

The exact solution ( |32] i for the generating function can in principle be expanded in x and y to yield the probability 
distribution P,„ „ for arbitrary m, n. For instance, the system is empty with probability 

exp {ye^~^') 



0,0 



Ho 



1+7 e 



-1 ^-r' . 



(33) 



pey + £(0, f) 

with p - yl{l + y), see Eq. ( [TT| . The expressions for the clone size distribution are simple when n -Q, that is for the 
clones without post-mitotic cells. Expanding F{x, 0, f) in powers of x and using F{x, 0, f) - 2m>o ^Pm,oit) we obtain 

(34) 



pexp(ye^' + f) „ 
"mfl = — — : — — — — P 



£(0, f) 



[pey + G(0,t)Y 

The probabilities P,„ „ quickly become very unwieldy for n > 0. 
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Figure 1: The permissible range in the model parameter space is < f < 1/2 and < y < oo. Shown are the curves [r= 1/4, r = 7/(1+ y)2 thick 
blue, and y = 2r/( l - 2r) thin green] where the solution has somewhat simpler forms. At the intersection of these three curves (y = 1, r = 1/4, □) 
the solution is particularly neat (see section|5j. The experimentally measured parameter values (y = 0.28, r = 0.08) for mice epidermis in jClayton| 
|et al.|[2007^ are depicted by the • symbol. 

5. Explicit results at the special point 

At the special point jc = l,rc - 1/4 in the parameter space we can solve everything explicitly. Indeed, in this case 
-y = and ( [32] ) becomes 

1 

Fix,y,t) = 1 + {l-y)e'' + 



(x + y-2)-^ -t/4 

Let us first extract the reduced distribution. Writing x - z,y = z and noting that 



j>0 



(35) 



(36) 



we conclude that 



G{z,t) = 1 +(1 -Z)e-'- y + y 



Expanding the latter expression in z around z = we get 



t t 

1 + ?r - 
2 2 



(37) 



In the scaling region 

equation P8c|l acquires a scaling form 



Ho = 1 + e"' - 



f + 2 



(f + 2)2 



it + 2)2 



f + 2 



s>2 



s ^ oo, t ^ oo, - — finite 
t 



(38a) 
(38b) 

(38c) 

(39) 
(40) 
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Recall that the exact expression (|6]l for the distribution of A cells also acquires an asymptotic scaling form; in the 
present case r - Vc - 1/4 it is given by 



Generally, by expanding 



and, for (w, n) (0, 0), (0, 1), 



P,n{t) 

we obtain 



16 



(f + 4)2 



f + 4 



m-1 



16 



-Imjt 



Po.i = -e-' + 



(f + 2)2 



p - 



(f + 2)2 



m + n 
m 



(41) 



(42) 



(43) 



2(f + 2) 

The probability that the system is empty is Po,o = Ho, so it is given by Eq. (|38a|. 

The clone size distribution greatly simplifies at this special point due to a mapping of our two-type branching 
process onto a single-type critical branching process. Indeed, at y = l,r - 1 1 A, the process can be reformulated as 



C 
C 



CC 




at rate 1/2 
at rate 1/2 



(44) 



where we assign the type A or B to each cell independently with probability 1 /2. This mapping holds if also initially 
we have an A or a B cell equiprobably. If the initial cell is type A, then from the solution for a single initial B cell ( [T7] i, 
and from the solution (|6]) of ( |44| i, we recover the behavior (|42]i-(|43]l due to the linearity of the problem. 



6. General results 



In sectionBlwe have found an explicit, exact expression for the generating function, equation ( |32| i, which is valid 
on the curve p3] l. This curve misses the parameter values (y - 0.28, r - 0.08) experimentally measured in mice 
tail epidermis ( Clayton et al. 2007 1, see Figure [T] In different tissues the parameters will probably take different 
values, so it is desirable to possess a solution in the whole range of parameters, i.e. in the strip < r < 1/2 and 
< y < oo. Surprisingly, using the reduction of the Riccati equation to the Sturm-Liouville equation we can find a 
general solution. 

We start with the general backward equation ([T9]l which we re-write in a canonical form 

F' = AF^ +BF + C 



The coefficients of the quadratic polynomial on the right-hand side are 



r 

yu ' 



B 



1 



1 +2r- 



C 



yu 



-d-uf 



(45) 



(46) 



To transform the Riccati equation ( |45] l into the Sturm-Liouville equation we perform the standard procedure ( |Bender| 
|and Orszag]|1978| l, namely we write F{u) as 



F = -— = (l°g^)' 
Az A 



(47) 



After this transformation, the first order nonlinear equation ( |45| ) turns into a second order linear differential equation 

z" + az' +liz^Q (48) 

with 



a - — 



^' I ^\ y-«-2r(l-M) 



uy 



P^AC^ 



r(l -u) 
uy 



(49) 



Now in ( [48| l the first derivative can be cancelled by writing z - OZ, with the condition O' — which leads to 

d) ^ g-5/"a(«'y«' (50) 

Then ( |48] l becomes a Shrodinger equation for Z(m) 

/4r- 1 y(l-2r)-2r 1 \ 
Z" + — ^ + ^ , + — rZ = (51) 

Equation ( |5T| resembles the Whittaker equation. Re-scaling the variable u and making changes in notations 

g = — , V = VI - 4r, w - (52) 

y 2yv 

we recast equation ( |5T| into a canonical Whittaker differential equation 

''^W4---A)z = (53) 



I 4 g 4^2, 
Its solution, up to an irrelevant constant factor, is 

Z(8) = M„,,o(g) + CW,,,o(g) (54) 

where M and W are the Whittaker functions (jGradshteyn and Ryzhik 2007^ , and C is a constant to be determined 
from the boundary conditions. 

Now we have to re-express the solution of Eq. ( |54] i in terms of the original variables. Following the steps that have 
been made, but backwards, we obtain 

F(«).-M^.Z^[logZ(«).log<l>(«)]' (55) 
A r 

Using [logO(M)]' = -all, see Eq. dSOl, we get 



'^^"^^7-M.,o(,).CW.,o(,)"^-""^ ^''^ 



Noting that 

(g - 2w)M„,o(g) -H (1 -H 2w)Mi+„,o(^) 



K.,o(^) = 



2^ 



(57) 



2^ 

we simplify ( |56] l and arrive at our main result 

m(1 +v)-r(l +2w) 

F = 1 - M H 

2r 

(1 + 2w)Mn.H,,o(g) - 2CWi-,v,.o(g) 
^ 2r ' M„,o(g) + CW„fl{g) 



(58) 



Recall that the parameters g, v, w are given by ( [52| l, and m = ( 1 - 3')e The constant C in Eq. ( [58] l is determined from 
the initial condition, F(m = 1 - 37) = x, to give 

^ ^ -gM,,o(g) + (l +2w)Mi+.,,o(g) 
0W„,o(g) + 2Wi+vv,o(g) 
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Figure 2: Exact survival probability S(t) = 1 - F(0, 0, given by {61) , as a function of time at the experimentally measured parameter values 
y = 0.28, r = 0.08. The dashed line is the large time asymptotic 1/rt. 



Here we introduced two more shorthand notations 



1 + 2w - £ + 



2r(x — y) + y — 1 

r 



(l-y)v 

r 



(60) 



The distribution of the total number of cells lls(t) can be obtained from G{z, t) - F(z, z, f). The survival probabihty 
of the cells at time t is 

5(f) = 1 -mo, f) (61) 



where F is given by ( pS) and ( [59) . Note that in computing F{x - 0,y — Qj) all the parameters that contain x and y 
simplify. Setting x = y = we get u - e^^' , 9 - 1 + 2w - (v + l)/y, and g - vjy. 

From the generating function ( |58| ) one can easily extract the clone size distribution Pnj„(t) or Tlsit) numerically. 
Let us start with the simpler total cell distribution Tlsit). The probability n5(f) is the coefficient of the z" term in the 
power series of G(z, t) as given by (l36|. One way to extract n5(f) is by using Cauchy's integral formula 



(62) 



Here the contour C goes counterclockwise around the origin in the complex z plane, within the radius of convergence 
of G(z) (we omitted the time argument for brevity). Consider a contour of a circle of radius R, and divide the circle 
into equal parts. Now the above integral (|62]i can be approximated as a sum 



^ k=0 



sIti/N 



(63) 



which is the discrete Fourier transform scaled by /? /N. This transformation can be performed incredibly efficiently 
by the fast Fourier transform (FFT) method, which is implemented in most mathematical software. This method is 
discussed and error terms are approximated in ( Cavers| |1978[ ). Some care is needed to choose the value of R to avoid 
numerical problems, as discussed in ( Cavers 19781). In our case the choice R - I was sufficient in all examples we 
considered. We can check the quality of this numerical method at the special point y = 1, = 1/4, where the explicit 
solution for P,„,«(0 is known ( |43] l. For example at f = 1, with N - 32 the numerical result for Ils(t) differs less than 
lO"'*" from the exact expression for s < 31, and it is precise to at least ten digits for s < 15. 
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For the full distribution Pm,n(t) one needs two separate contour integrals in both the x and the y planes, which then 
leads to applying the discrete Fourier transform A^^ times. We have checked the results against the numerical solution 
of the forward equations (|2]) and found a perfect agreement (up to about 7 digits). This method has been used to obtain 
our figures [3] and |4] for Tlsit) and Pm,ni^)^ respectively. In (Clayton et al. 20071 the initial cell is considered to be A 
or B with certain probabilities. The corresponding probability distribution is then a simple linear combination of the 
distribution Pm^nif) we just obtained and the trivial distribution resulting from a single initial B cell ( fTTj i. Since B cells 
just die at a fixed rate, their only effect (apart from f 0,0 and P\fi) is to rescale Pm,n{t)- 

Thus we have obtained exact results (|58]l-(|59]l for the generating function, which can be easily transformed back 
to probabilities. Moreover, these exact results simplify in a few special cases (Appendix [A]i and in the scaling limit 
(Section|7]i. 



7. Scaling limit 



In the large time limit, the distributions Pn,m(t) and Ils(t) simplify. Let us consider first the reduced distribution 
G(z, f). In the large time limit the interesting range of s is s ~ f, see e.g. ( (39| ), and therefore the interesting range of z is 
{ \ - z) ~ ~ . Hence we consider the r — > oo,^ — > 1 limit with ^ - rt(\-z)lp kept constant, where p - y/(l +y), 
see ( [TT] l. 

In order to perform the scaling limit we need the following small argument (jc <K 1) limits of the Whittaker 
functions 



(64) 



r(i/2- w) 

Here F is the gamma function, ^{z) = F'(z)/F(z) is the digamma function, and -y^ - 0.5772 ... is the Euler constant 
(Gradshteyn and Ryzhik 2007 1. We also need the identity for the digamma function (Bender and Orszag 1978 1 



1 



^ _ _ w -I/. - - w = 



1 



1 +2w 



(65) 



Taking the z ^ 1 limit of the constant term (|59| is particularly easy, since it is independent of time. We find 



C = Cod - z) + 0(1 - z)''\ Co = r 



I + y 



r(l/2- w) 



(66) 



Now we substitute this expression into F(z, z, t) of ( |58| ), using 



y r(l + y)t 
Some care is needed with terms of type Wwfiig), where in 

logg = -rf + log 



r(l + y)t 

a term proportional to t appears. In the first order of 1/f we obtain 

1 ^ 



G(^,0 = 1 



rt (+1 



(67) 



(68) 



(69) 



In the s, f — > 00 scaling limit with constant yu - ps/rt , the generating function G(^, f) of ( [36| l becomes a Laplace 
transform of 11 



G(^,f)^ - n,(^)(f)e-*fl'/z 
P Jo 



(70) 
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hence we can obtain the asymptotic limit of the probability nj(f) by an inverse Laplace transform ( Gradshteyn and 



Ryzhik 2007 1 



-n,(f) ^ r'iGiCM = -e-^' + (l - -\6{ti) (71) 
p rt \ rtj 

The first term describes the distribution of the surviving cells, while the second term stands for the extinction of cells. 
Consequently, the large time survival probability of the population is 1 jrt, or the extinction probability no(f) ~ 1 - 1 / rt. 
In Figure |2] we plotted the exact survival probability Sit) — 1 - no(f) of ( |6T] i together with the large time asymptotic 
Ijrt. The above asymptotic results of course agree with the explicit results in the special point of Section|5] 
From ( |7T| , the regular part of the distribution HAf) can be written in a scaling form as 

n.(0 = 4exp(-^)=^na.) (72) 

with the time independent scaling function 

nOu) = e-^' (73) 

We demonstrated this scaling in Figure |3] where the exact expressions ( (58] l for 115(0 depicted for different times 
as a function of the scaling variable ^ - psjrt, and the values converge to the scaling limit ( |73| ). Note that this scaling 
limit has been already guessed in ( [Clayton et al.|[2007] l, and derived in ( [Klein et aL] [2007 ) in the realm of continuous 
approximation, that additionally assumed that the B cell population remains "slave" to the A cell population. We will 
see that the latter, potentially uncontrolled approximation is not merely appealing, it is asymptotically correct. This 
will become evident from the full distribution. 

Similarly to the total cell distribution, we can also obtain the scaling limit of the whole Pm,n(f) distribution from 
( |58| ). Taking the r — > 00 limit while keeping ^ = n(\ - x)lp and 77 = rt{\ -y)/p finite, up to first order in 1/f we obtain 

1 £y + ri 

F(^,rj,t)^l--- (74) 

Of course F{(, t) - G{(, t) of (|69]l in this limit as well. Now we need to perform a double inverse Laplace transform 
to obtain Pm_„(f) as a function of mjt and n/f, in the limit m, n,t ^ 00. The extinction probabiUty is again Po,oit) - 

11 



> 



0.025 



0.02 - 



0.015 



0.01 



0.005 - 




(m-Yn)(p/cort) 



Figure 4: The exact probability Pm,n(t) of having m type A, and n type B cells at time t, as given by \58\ . This probability is depicted at different 
times in terms of the scaling variable y = (m - yn) ^p/airt, at /j = (m + n)p/rt = p/r a: 2.7. The points collapse on the Gaussian limit curve jSTj. 



no(f) - 1 /rf in the first order of 1/f. The probability Pm,n{t) for m,n > in the scaling limit becomes 



2,3 S'^P 



p(m + n) 



rt 



(75) 



Hence in this limit there are precisely y times as many A cells as B cells, while the distribution of the cells is given by 
( |72l ). According to the experiments of |Clayton et al. ( 2007 1 in skin tissue, where y = 0.28, the model predicts about 
four times more post mitotic cells than progenitor cells in a clone for large times. 

It is possible to give a more detailed description of the cell distribution by taking a different large time limit, 
namely we need to take the limit n, m, f ^ oo in such a way that the following fractions are finite 



= 0(1), -=0(1), 



m — yn 



-0(1) 



(76) 



This limit reveals the "shape" of the Dirac delta in ( |75] l. This asymptotic limit is of course encoded in the exact 
results for the generating function ( |58| l. Unfortunately, to extract the asymptotic is far from straightforward. Indeed, 
even from a simple expression for the multivariate generating function, it is usually extremely difficult to extract 
the asymptotic of the coefficients (let alone the exact expressions for the coefficients). This situation is perhaps 
surprising as in the univariate case there are various techniques, the most powerful is the use of complex analysis 
and the saddle point method. In the multivariate case, the usage of complex methods is much more limited and 
challenging; for recent progress, see ( Flajolet and Sedgewickj 20091 and (Pemantle and Wilson 2008| l. In our case, 
there is an additional difficulty as the explicit expression for the generating function is not a simple rational function 
as e.g. in most examples in (jPemantle and Wilson 2008), but it involves the Whittaker functions. Hence instead of 
extracting the scaling limit from the exact solution, we outline another approach in the next section that also shows an 
independent way of handling the problem. 



8. Fokker-Planck approximation 



We shall use a more direct procedure which is however approximate, for instance it does not even provide the 
asymptotically exact value, 1 / rt, that the clone size is non-zero. However, up to this amplitude one can obtain an 

12 



expression for the probability distribution Pm,«(0 which is typically asymptotically exact in the scaling region ( |76| l. 
The method is essentially the Fokker-Planck or diffusion approximation (? van Kampenj |1997| l. One starts with the 
master equation (|2|, and treats m, n as continuous variables. This should be valid when m, n » 1. In this region one 
can further expand the right-hand side of Q in the Taylor series to give (we shortly write P instead of Pm,n) 

d 1 

(m - l)Pm-i,n = mP- —mP+ - —r mP + . . . 

dm 2 dm^ 



d 1 d- 

mP,„^„-i = mP - -^mP + -— mP + 
on 2 on-^ 



(77) 



d 5 1 

(m + l)P,„+i n-2 - mP + — mP - 2— mP H mP - 2- — — mP + 2 — - mP + ... 

am on 2 om'^ oman on-^ 

d 1 52 

(n + l)P,„,„+i ^nP+ —nP + -—rnP + ... 

on 2 on^ 

Using these expansions and ignoring the higher order terms we turn the master equation into a partial differential 
equation 



dP , ^ dP ^ dP 

— - yP + (y - 2r — m + yn) — — h 2r - — 
ot on dm 

_^ d^P 2 ^ (1 + 20?" + rn d^P 

dm^ dmdn 2 dn^ 



(78) 



which is the Fokker-Planck equation in our problem. 
Let us change m, n to the variables 

s = m + n, 6 — m — yn (79) 

The Fokker-Planck equation becomes 

dP dP dP 

— ^yP + {y-S)— + [2r + y{6 + 2r- y)] — 
at as do 

Q2p Qlp Qlp 



with 



^^ 2ysUl-y)6 ^ 
2(1+7) 

2(l+2r)(ys + 6) + y(s-6) ^.7S + 6 

C - y — h r( 1 + 2y) ■ 



(80) 



2(1 -H 7) 1+7 
Since 5 <s; s in the scaling region (|76]l, the above coefficients simplify to 

A = pi, B = yps, C =ps [r(l + yf + 7^] (81) 

We akeady know the dependence on s, namely P ~ exp (-77)- To determine the dependence on 6 we keep only the 
dominant terms in the Fokker-Planck equation (l80| and obtain 



dP d^P 



Note that all terms in (l82| are of the order of P: 



dP P d^P P 
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The latter estimate follows from s ~ 6^ and it actually explains the choice 6 ~ f'^^ in the scaling region ([76|. Note 
also that the neglected terms from the Fokker-Planck equation ((80| are indeed sub-dominant, e.g. 



glp p p d^P P P 

A—r~s^^-, B —— ~ s — ^ - (84) 

Solving ( [82) l, which is essentially an ordinary differential equation with respect to 6, we find 

P ~ exp (-— ) , with w = 2(r + p^)(l + y) (85) 



(OS 



Therefore the full scaling solution reads 



P„,M = . exp I - — I (86) 



{rty s/jmjs \ rt ojs 

The amplitude, including the 1/ factor, is obtained by requiring I!., = J 6)d6/(l + y), using ( |72| l. Note that the 
distribution ( |86l ) is normalized as J Pdmdn - J P ds d6 /(I + y) - (rf)"'. 
The limit distribution ( |86l ) can be written in a scaling form 

Pm,nit)^ ri^J-P(M,yl with P(f^,v)^^— (87) 

(r?) ' V w 

fi - — , V = (88) 

rt V ^^ff 

This scaling is probed in Figure |4] using exact values for Pmj,{t) from ( |58] l. One can see that the scaling limit ([86| 
provides an excellent approximation already for times f > 10, and the finite time curves converge to the scaling 
function ( (87] l. Note also that in the special point y - 1, r = 1/4 the distribution (|43j converges exactly to the scaling 
limit dSTli. 



with scaling variables 



9. Discussion 



We derived an exact solution for a two-type branching process. We investigated a specific stochastic process that 
has been proposed to describe measurements of murine tail epidermis ( Clayton et al.|[2007] l. The chief ingredient 
of the stochastic process ([TJ is the self-duplication and differentiation of the progenitor cells without measurable 
contribution from stem cells. (Stem cells activate during repair from severe injuries.) The same mechanism apparently 
underlies the maintenance of pancreatic islets ( Dor et al. 2004) and lung homeostasis (Giangreco et al. 2009!). 

An exact solution of the specific two-type branching process ([T} raises the hope that other two-type branching 
processes could be amenable to analytical treatments. Some two-type branching processes have been suggested 
long ago in the context of tumor formation (Kendall 1960). Indeed, cancer is often arises when a progenitor cell 
undergoes a series of mutations in a way that the proliferation of a mutant clone dominates the differentiation or death 



(Fearonetal. 1987 Fearon and Vogelstein 1990 Dingli et al. 2007a Nowak 2006 


Attolini and Michor 


2009|l. 


The complication is that cancer typically involves multiple mutations ( 


Armitage and Do 


1 1954 Beerenwinke 


I et al. 



The prominent feature of our analysis is the disregard of spatial characteristics. In the context of epidermis, one 
might want to consider the two-dimensional version of the two-type branching process ([T]). The spatial model is partly 
amenable to analysi s (Klein et al.^^2008) due to an intimate connecti on with models of voting and monomer-monomer 



catalytic reactions (Krapivsky 1992 Frachebourg and Krapivsky 1996[ Liggett 19991. Intriguingly, although the 



model presented in this paper completely disregards real space, it already provides excellent fit to experimental data. 
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A. Special Cases 

As a check of self-consistency it is useful to extract the explicit results of section [5] from the general approach 
of section [5] At the special point y= 1,'"= 1/4 the potential in equation ( |5T| ) is purely quadratic V = 1/4m^, and 
the solution of ( |5T| becomes Z = ^[u{C + log u). After transforming Z(m) back to F{t) and fitting to the boundary 
conditions, we indeed recover ( |35| l. 

Exact results (|58])-((59| also simplify on a few lines in the parameter space. 



A.l. Horizontal line r — 1/4 

On this line the M-independent term in the potential in Eq. ( |5T] l vanishes and the Shrodinger equation becomes 



z".i(z^«->.«-)z.o 



(89) 



In the M ^ limit, Z(m) behaves as ■\/u. This suggests to choose k - y/U as the basic variable and seek solution 
proportional to k. Hence we write 

Z(m) = kG(k), u = ^— k^ (90) 
1 -7 

The amplitude y^/(l - y) has been chosen to get rid off y in the coefficients of the governing equation for G(fe): 

(fG ^ \ dG ^ _ ^ ^^j^ 
dk^ k dk 

Solutions to this equation are linear combination of the modified Bessel function I^ik) and Koik), i.e. 

G(k) = CMk) + C2Ko(k) (92) 

Then the function F{u) is given by 

CIi{k)-Kiik) 

F^\+u + 2yk )' ' (93) 
ClQ{k) + Ko(k) 

where we used identities I'^ix) = Ii(x),K'^{x) - -Ki(x) and C = C1/C2. We can re-write this as 

CIi{k)-Ki(k) 

Fix, - 1 + (1 - y)e-^' + 2yk ^ ' , ^ ' , (94) 

Ck){k) + Ko(k) 

with 

k = y-'e-y"^ Vd-rXi-y) (95) 

where C = C(x, y) is determined by matching to the initial condition F{x, y,t -0) - x. One gets 

yKKd'<)-(l-^^)KoiK) 
C = (96) 

r^/i(^) + (i-:^)/o(^) 

with 



^ = ^(f = 0) = y-' V(i-r)(i-y) (97) 

A.2. Special curve y — 2r I {\ — 2r) 

On this curve w -Q, hence the term proportional to g ' vanishes in the Schrodinger equation ( |53] l, which then can 
be solved in terms of Bessel functions. Here instead, we derive the simplified form from the general solution ( |53| l. 
For this we need some limit properties ( [Gradshteyn and Ryzhik|[2007| l of the Whittaker functions 

Mods) = <ghigl2) 
Wo,o(g) = ^/^Ko(g/2) 

Muo(g) ^(1-8) yl8lo(8/2) + g''^Ii(8/2) (98) 

WM = [-(1 - ^) y/gKo(g/2) + g''^Ki(8/2)] 
15 



By using these expressions in ( |58| ), we obtain 




h{gl2)-CK,(gl2) 
UgH) + CKo(g/2) 



(99) 



with constant 



C = 



-Xk{g/2) + vh{g/2) 
xKo(g/2) + vKi(g/2y 



X = 



2r{x — y) 



- 1 



(100) 
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